Theoretical analysis of dynamic chemical imaging with lasers using high-order 

harmonic generation 
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We report theoretical investigations of the tomographic procedure suggested by Itatani et al. 
[Nature, 432 867 (2004)] for reconstructing highest occupied molecular orbitals (HOMO) using high- 
order harmonic generation (HHG). Using the limited range of harmonics from the plateau region, 
we found that under the most favorable assumptions, it is still very difficult to obtain accurate 
HOMO wavefunction, but the symmetry of the HOMO and the internuclear separation between the 
atoms can be accurately extracted, especially when lasers of longer wavelengths are used to generate 
the HHG. We also considered the possible removal or relaxation of the approximations used in the 
tomographic method in actual applications. We suggest that for chemical imaging, in the future it 
is better to use an iterative method to locate the positions of atoms in the molecule such that the 
resulting HHG best hts the macroscopic HHG data, rather than by the tomographic method. 

PACS numbers: 42.65.Ky, 33.80.Rv 



I. INTRODUCTION 

In the microscopic world, chemical reactions and bi- 
ological transformations occur on the time scale of the 
order of picoseconds or less. Understanding the inter- 
mediate steps of these reactions or transformations has 
been the dream of physical and chemical scientists [j], 0] • 
While X-ray and electron diffraction have served well to 
image the structure of large molecules, current technol- 
ogy limits such tools to time resolution to tens of picosec- 
ond to sub-picoseconds [!, H[ , making them less suitable 
for time-dependent studies. On the other hand, infrared 
lasers of durations of tens to sub-ten femtoseconds are 
becoming widely available today, thus it is natural to 
ask whether femtosecond lasers can be used for dynamic 
chemical imaging of molecules. Since the wavelength of 
infrared lasers is much larger than the interatomic spac- 
ings, these lasers cannot be used for diffraction measure- 
ment directly. Nevertheless, recent experiments Q have 
shown that the shape of the highest occupied molecu- 
lar orbitals (HOMO) can be revealed by measuring the 
alignment-dependence of the ionization rates by sub-ten 
femtosecond lasers. Similarly, limited experimental evi- 
dences and theoretical calculations so far have also shown 
that the high-order harmonic generation (HHG) can re- 
veal the structure of molecules IE 0, [1] , including even 
two-center interference 0, [13, lilf - Self-diffraction from 
the rescattering electron following laser ionization has 
also been suggested for imaging molecules [HI, [l]| flij . 
Despite many such discussions in the literature in the 
past decade, few real demonstrations have been achieved 
so far. 



they have reconstructed the highest occupied molecular 
orbital of N2 molecules accurately using the tomographic 
procedure from the measured HHG spectra. Using a 
weak laser beam first to partially align N2 molecules, 
a more intense short probe pulse (30 fs) was then used 
to generate high-order harmonics. The polarization di- 
rection of the probe pulse was varied from parallel to 
perpendicular direction with respect to the pump pulse. 
The retrieved HOMO wavefunction has been shown to 
be in surprisingly good agreement with the HOMO cal- 
culated from the quantum chemistry code. In a more 
recent paper, taking that the HHG are generated from 
a many-electron N2 molecule, it was further shown [l6j 
that the reconstructed wave-function records the image 
of the Dyson orbital plus exchange contributions from 
the inner shells. 

The success of Itatani et al. has generated a lot of 
interest. If it is generally applicable, it would offer the 
means for time-resolved dynamic chemical imaging at the 
resolution of ten to sub-ten femtoseconds. However, a 
number of skepticism have been raised recently flo ] . 
Rigorously speaking, in quantum mechanics, the ground 
state wavefunction is not measurable alone since any 
measurement would invoke final states. The tomographic 
procedure used by Itatani et al. relies on the approxima- 
tion that the continuum states of the molecule in the 
recombination process be represented by plane waves, an 
approximation that is well known to be invalid in pho- 
tochemistry, unless the continuum electron's energy is 
in the keV range. This is not the case for HHG in the 
plateau region where the photons are from a few tens to 
about 100 eV. 



In a recent paper, Itatani et al. 15[ reported that The starting point of Itatani et al. is the three-step 



2 



model of high-order harmonic generation [lj| [2(J, [HI- 
The harmonic yield is assumed to be given by the product 
of three terms, one is the initial tunneling ionization rate, 
the next is the propagation of the ionized electron, and 
the third is the recombination of the electron with the ion 
to emit high-order harmonics. By assuming that the first 
two steps are nearly independent of whether the target 
is molecular or atomic so long they have nearly identical 
ionization potentials, Itatani et al. were able to extract 
the recombination dipole matrix element for individual 
N2 molecules from the measured HHG in the plateau re- 
gion, as well as their variation with angles. These dipole 
matrix elements arc then used to obtain the molecular 
wavefunction using the tomographic imaging method. 

Two important approximations are used in this tomo- 
graphic procedure. First, the dipole moment of individ- 
ual molecules in the direction parallel and perpendicular 
to laser's polarization can be obtained from the measured 
macroscopic HHG. For this, they assumed that the har- 
monics in the direction parallel to the laser polarization 
satisfies complete phase matching condition. Experimen- 
tally this is done by placing the gas jet after the focus 
of the laser, thus eliminating the contribution to HHG 
from the long trajectories [23, HH. They did not mea- 
sure HHG in the direction perpendicular to the laser's 
polarization, but the dipole moment in the perpendicu- 
lar direction was assumed. In practice, if the molecules 
are fully aligned, i.e., all pointing in the same direction in 
space, then the measurement of macroscopic HHG in the 
perpendicular direction would allow the determination 
of the dipole moment in the perpendicular direction. For 
molecules that are only partially aligned, such as those 
aligned by a pump laser, it is not straightforward to ex- 
tract the perpendicular dipole moment of each molecule 
from the measured macroscopic HHG. 

The second important approximation in Itatani et al. 
is that the continuum wavefunctions in the recombination 
dipole matrix elements can be approximated by plane 
waves. This approximation is at the heart of the tomo- 
graphic procedure. Implicit to this also is that the dipole 
moment is available over the complete spectral range of 
the harmonics. This is not so since the useful harmonics 
are usually limited to the plateau region only. 

In this paper we examined the assumptions used in the 
tomographic procedure. We used the Lewenstein model 
[2TI ] to calculate the HHG's for fixed molecules with dif- 
ferent alignment angles with respect to the laser polar- 
ization. In this model, the dipole matrix element is cal- 
culated using plane waves for the continuum states. We 
thus take the calculated HHG from the Lewenstein model 
to be the "experimental" data, and follow the procedure 
of Itatani et al. to extract the "experimental" dipole ma- 
trix elements. If the three-step model is exact, the "ex- 
perimental" dipole matrix element thus obtained should 
be identical to the "exact" theoretical ones which are cal- 
culated using plane waves for the continuum states and 
the HOMO for the initial state. In reality, we found that 
the "experimental" dipole matrix element agrees with the 



theoretical one only in a limited spectral range - in the 
plateau region. With such a limited range, we found that 
the HOMO wavefunctions derived from the tomographic 
procedure to be approximate. Even though the accuracy 
can be improved by using HHG generated with lasers of 
longer wavelength to extend the plateau region, we were 
unable to obtain accurate HOMO wavefunction. 

In view of the failure in obtaining accurate HOMO 
wavefunctions, we next turn to the possibility of using 
HHG to extract some information on the structure of a 
molecule. We note that in chemical imaging it is not es- 
sential to reconstruct the wavefunction of a molecule. A 
great deal is understood if the positions of the atoms in 
the molecule can be determined at the intermediate time 
steps. We thus ask what partial information of a molecule 
can be revealed by the tomographic procedure under the 
most favorable conditions assumed in our model. Using 
lasers with wavelength of 1200 nm to generate plateau 
harmonics from N2 and O2 molecules at different inter- 
nuclear distances, we showed that the internuclear dis- 
tances obtained from the retrieved HOMO wavefunction 
is good to better than a few percents. Furthermore, the 
geometry and the symmetry of the HOMO, including the 
nodal surfaces, can also be retrieved. Such precision is 
more than adequate for dynamic chemical imaging and 
warrants further test for complex molecules where the 
positions of some atoms undergo large change, in pro- 
cesses such as isomerization or dissociation, and where 
the time scale of transformation is of the order of tens of 
femtoseconds. 

The rest of this paper is arranged as follows. In 
Section II we describe the theoretical methods and as- 
sumptions used for the present tomographic imaging 
of molecules using high-order harmonic generation of 
aligned molecules by lasers. Section III presents the re- 
sults of extracting HOMO wavefunctions and the inter- 
nuclear separations. In Section IV we revisit the assump- 
tions and address the intrinsic limitations of the tomo- 
graphic imaging method. We showed that these limita- 
tions are not easily removed in real applications, making 
it difficult to apply the tomographic method for chemical 
imaging in general. In Section V, we suggest alternative 
methods for imaging molecules with HHG, based on the 
iterative procedure. Assuming that the HHG for a sin- 
gle polyatomic molecule can be calculated accurately, the 
propagation effect and the effect of partial alignment of 
the molecules can be taken into account by the standard 
methods. Thus macroscopic HHG can be calculated for 
each configuration of atoms in the molecule which can 
then be compared to the experimentally measured HHG. 
Since HHG can be generated for different alignment an- 
gles of the molecules using lasers of different wavelengths, 
an optimum atomic configurations can be derived by the 
set that best fits the experimental HHG. We note that 
in typical chemical reactions, only a few atoms change 
positions significantly. Thus the number of parameters 
that have to be determined is small. We then give ex- 
amples of the isomerization of acetylene C2H2 and of hy- 
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drogen cyanide HCN to show the alignment dependence 
of the HHG of these molecules under isomerization. In 
the last Section we summarize our conclusion on the to- 
mographic procedure and discuss future plan for study- 
ing dynamic chemical imaging with lasers. Atomic units 
m = e = h = 1 are used throughout the paper, unless 
otherwise indicated. 



II. THEORETICAL METHODS 

A. The generation of HHG using Lewenstein model 

In this paper, in order to generate the HHG spectra 
from molecules with fixed alignments, we use the Lewen- 
stein model, as extended by Zhou et al. ((J. Without loss 
of generality, we assume that the molecules are aligned 
along a;-axis in a laser field E(t), linearly polarized on 
the x-y plane with an angle 9 with respect to the molec- 
ular axis. The parallel component of the induced dipole 
moment can be written in the form 

r°° / 7r \ 3 ^ 2 
D\\(t) = ij dT^-^-j^) [co S 0d*(t)+sin0^(t)] 

x [cos 9d x (t - t) + sin 9d y (t - r)]E(t - r) 

x exp[-iS st (t, r)]a* (t)a(t - r) + c.c. (1) 

where d(t) = d[p st {t,T) - A(t)], d(t -r) = d\p st {t,T) - 
A(t — t)] are the transition dipole moments between the 
ground state and the continuum state, and p s t(t,r) = 
J. A(t')dt' /t is the canonical momentum at the sta- 
tionary points, with A the vector potential. The per- 
pendicular component D±(t) is given by a similar for- 
mula with [cos9d x (t)+sm9d*(t)} replaced by [sin#d*(t)- 
cos 9d* y (t)} in Eq. (1). 

The quasiclassical action at the stationary points for 
the electron propagating in the laser field is 

where I p is the ionization potential of the molecule. In 
Eq. ([T]), a(t) is introduced to account for the ground state 
depletion. 

The ground state electronic wavefunctions were ob- 
tained from the GAMESS code J24J- Within the Sin- 
gle Active Electron (SAE) approximation, we take the 
HOMO for the "ground state" , the transition dipole d(k) 
is given as (k\r\^o(r)) . The continuum state is approx- 
imated here by a plane- wave |fe). This model has been 
shown to be able to interpret the experimental results 
so far for N 2 , 2 @, H and C0 2 [Tfl, HE S3- In this 
paper we neglect the effects of depletion, i.e., we assume, 
a(t) — a*(t) — 1. Laser peak intensity of 2 x 10 14 W/cm 2 
and duration of 30 fs (FWHM) are used throughout this 
paper. 



B. The procedure of tomographic imaging 

In the tomographic procedure, suggested by Itatani et 
al. [IB], the intensity of the emitted HHG is approxi- 
mated by 

S{uj, 9) = N 2 (9)oj 4 \a[k{uj)}d(uj, 9)\ 2 . (3) 

Here d(u>,0) = (*f>o(r; 9)\r\ exp[ik(u)x]) is the transition 
dipole between the HOMO and the continuum state, ap- 
proximated by a plane- wave with/c(w) = (2UJ) 1 / 2 -, a[k{uj)} 
is the amplitude of the continuum plane wave of returned 
electrons; and N{6) is the number of ions produced. 
^o(r;9) represents the HOMO, rotated by angle 9 (the 
Euler angle between molecular axis and the laser polar- 
ization direction). Note that in Eq. §5§ perfect phase- 
matching was assumed in the experimentally measured 
spectra. This equation can be understood as a further 
approximation to the three-step model [28| . In the case 
of strong-field tunneling ionization, it can be assumed 
that the ionized electron wave packet depends weakly 
on the molecular structure, but rather on the ionization 
potential only. Thus, one can get rid of a[/c(w)] by mea- 
suring the HHG from the reference atoms with the same 
I p (for N 2 the reference atom is Ar, for 2 it is Xe). 

From the measured HHG and ionization data one can 
extract the transition dipole d(k, 9) by using formula 

|d(M)l = N(9y 1 \d ref (k) | y/s(u>, 9)/S ref M (4) 

with d re f(k) being the transition dipole moment of the 
reference atom. From Eq. Q we can extract only the 
amplitude of the dipole moment. However, it is well- 
known [291 ] that the minimum of harmonic signal and a it 
phase jump occur at the same harmonic order. Following 
Itatani et al., we assume that the dipole changes sign 
when its absolute value goes through a minimum close 
to 0. 

In this study we will start with HHG from individual 
molecules which were calculated using the Lewenstein 
model. This allows us to bypass the issues related to 
phase matching and other macroscopic effects for the 
time being. Following Itatani et al. we obtain the d(u>, 9) 
from Eq. © except replacing N 2 {9) in the equation 
by N(9) (since phase matching is not used for a single 
molecule). The ionization rates N{&) are calculated using 
the molecular tunneling theory (MO-ADK) [33 ], 

One major difference from Itatani et al. is that in our 
calculations we use the dispersion relation 

k = ^2{lo-I p ) (5) 

with I p being the ionization potential. The question 
about including or not including I p in Eq. (JSJ) was dis- 
cussed in a number of papers (see, for example, [30l l3ll] ) : 
it is a question of how to relate the energy of the re- 
colliding electron with the energy of the emitted photon. 
For radiative recombination or photoionization, Eq. ([S]) 



4 



is the correct relation with k being the momentum of 
the electron in the asymptotic region. In the Lewenstein 
model the re-colliding electron does not see any potential 
from the ion so k 2 /2 is also the kinetic energy of the elec- 
tron near the nucleus. Indeed we found that the deduced 
dipole moment from HHG calculated using Lewenstein 
model works well only when Eq. ^ was used. In Itatani 
et al, the relation k = \/2cu was used, but their dipole 
moment was extracted from the experiment. 

From the dipole moment parallel and perpendicular to 
the polarization direction, we first rotate the two compo- 
nents to the molecular frame. The wave function is then 
obtained by using the Fourier slice theorem [33[. The 
basic formulae for tomographic procedure 

xW(x,y) = /* d9 J + °° d^fcC* cos 0+2/ sin0) 

x [cos 9d x (cj, 0) + sin 9dy (to, 9)] (6) 

yV(x,y) = d6 f+°° duje lk< - xcose+ y sin ^ 

x [- sin 0d x (ui, 9) + cos 9d y (uj,9)] (7) 

can then be used to obtain the HOMO in the molecular 
frame. Note, here and in the following for the wave- 
function integrated along the third direction we use the 
notation 

/+oo 
*o(x,y,z)dz. (8) 
-oo 

In principle, the two equations ^ and ([7]) are equiv- 
alent. In practice, however, incomplete dipole moment 
data would give different results from equations ^ and 
0. We take the results to be the mean from the two 
equations. If we do not know the symmetry of the 
molecule the HHG should be measured over the whole 
alignment range, from 0° to 180°. For homonuclear di- 
atomic molecules the symmetry allows us to consider 
data from 0° to 90° only. Furthermore, from the in- 
vestigation of tomographic transformations © and ([7]) 
using "exact" theoretical dipole moments we found that 
an angular step of A9 = 5° is good enough for extracting 
the wave function accurately. Thus, instead of integra- 
tion over the whole range we use the summation over the 
19 angles from 0° to 90° only. 

III. RESULTS AND DISCUSSION 

A. Alignment dependence of HHG for N2 and O2 

In this paper, we generate HHG spectra using the 
Lewenstein model, extended for molecules, as described 
in Sec. H(A). In Fig. 1, we show the HHG intensities 
along the direction parallel to the laser polarization for 
N2. Here the molecules are assumed to be fixed in space 
and the directions of laser polarization are varied from 
to 90° (in step of 10°) with respect to the molecular 




HHG order 

FIG. 1: (Color online) HHG spectra with parallel polariza- 
tion (with respect to the laser polarization direction) from N2 
aligned at 10 different angles between and 90°. The laser 
wavelength of 800nm and intensity of 2 x 10 14 W/cm 2 are 
used. 



axis. The laser wavelength is 800 nm. Note that our 
results are for single molecules and only odd harmonics 
are used. For each (2n + l)-th harmonic, we integrated 
the calculated harmonic spectrum from the (2n)-th to 
(2n + 2)-th harmonics. For the tomographic imaging we 
need the harmonic yields in the direction perpendicular 
to the laser polarization as well. Figures 2 and 3 show the 
alignment dependence of some selective harmonic orders 
for N2 and 02- Both polarizations are presented. For 
parallel polarization, the HHG signals are peaked near 
parallel alignments for N2 and about 45° for O2. For 
perpendicular polarization, the HHG signals of N2 are 
largest near 45°. 

To compare the HHG calculated here directly with 
those measured by Itatani et al. [l5| is difficult. Their 
data are from the macroscopic medium such that prop- 
agation effect is included. Furthermore, the molecules 
from the pump beam are only partially aligned. Nev- 
ertheless by assuming phase matching and that the 
molecules are fully aligned, they extracted the dipole mo- 
ment of the individual molecules from the macroscopic 
HHG. In Itatani et al. the HHG in the perpendicular- 
direction was not measured, but maximal HHG inten- 
sity near 45° was used implicitly [H]]. This assumption 
is consistent with our calculated results for N2 shown in 
Fig. 2. For O2, our calculations show that the perpen- 
dicular HHG intensities almost vanish near 45° and peak 
at two intermediate angles near 30° and 70°. 

In reconstructing the HOMO's, we typically use 19 
alignment angles from to 90° (in step of 5°). This is 
also the typical setup used in the experiments by Itatani 
et al [HI]. We found that using a finer angular grid do 
not improve the results. 

As we will show, in order to improve the quality of the 
reconstructed HOMOs, one would like to use the HHG 
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Alignment angle 9 (degrees) 

FIG. 2: (Color online) Alignment dependence of selective 
HHG orders for parallel (top) and perpendicular (bottom) 
polarizations from Na. The laser parameters are the same as 
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FIG. 3: (Color online) Same as in Fig. 2, but for O2. 



over an extended plateau region. It is better not to use 
higher laser intensities, since it would result in depletion 
of molecules to distort the simple model [Eq. (3)]. In- 
stead, we choose to use a laser of longer wavelength, say, 
of 1200 nm. We chose the same peak intensity of 2 x 10 14 
W/cm 2 , as the 800 nm laser. In Fig. 4, we present HHG 
spectra from N2 aligned at 45° , for 800 nm (black curves) 
and 1200 nm (red curves). Clearly, the plateaus are ex- 
tended to much higher harmonics. We note, however, 
that one harmonic order is 1.55 eV for 800 nm laser and 
1.03 eV for 1200 nm laser. 




Harmonic order 

FIG. 4: (Color online) HHG spectra from N 2 aligned at 45° 
for 800 nm and 1200 nm lasers. The laser intensities of 2 x 10 14 
W/cm 2 are used in both cases. 



B. Extracted dipole matrix elements 

Firstly we use the HHG spectra obtained in Sec. Ill (A) 
(to be called "experimental" in the following) for N2 in 
the case of 800 nm laser to extract the dipole moment 
in the laboratory frame by using Eq. Ar(3p) (with 
Ip = 15.76 eV) is used as the reference atom, as sug- 
gested by Itatani et al [Hj]. The data are available for 19 
angles 9 in the range of [0 — 90°] . These extracted "ex- 
perimental" dipole moments will be compared to those 
calculated theoretically. For demonstration, in Fig. 02a) 
we show the "experimental" and theoretical results of 
d x (k,6) for the angles 9 of 0°,30°,60° and 90°. Simi- 
lar results are shown in Fig. HJa) for the y-component of 
dipole d y (k, 9). In order to compare with the theoretical 
data, we have normalized the x-component of the dipole 
at k 2 — 2.5 for 9 = 0, The normalization factor is then 
fixed and used for all the other alignment angles as well 
as for y-components. This normalization factor is also 
used for N2 with the 1200 nm laser. 

From the Figs.[5]Ja) and EJa) we see that the "exper- 
imental" dipoles compare well with the theoretical ones 
in the range of k 2 = [0.9 — 4.2], which is equivalent to 
the useful range of HHG from H17 to H47 for the 800 nm 
laser. The lower limit H17 is due to the inaccuracies in 
the Lewenstein model for low harmonics, while the upper 
limit H47 is due to the harmonic cut-off. For the 800 nm 
laser and intensity of 2 x 10 14 W/cm 2 the cut-off position 
is at about H35 for N2 , as can also be seen in Fig. [TJ 

We can enlarge the useful region of HHG by using 
longer wavelength lasers. In Figs.[5fb) and[^b) we show 
the "experimental" dipole moments in the case of 1200 
nm laser with the same intensity. The useful dipole mo- 
ments now cover the range k 2 = [1.0 — 7.5], or H27 to 
H113, of the 1200 nm laser. We note that the agree- 
ment between the "experimental" and theoretical dipoles 
appears to be better for higher plateaus. This is not en- 
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FIG. 5: (Color online) x-component of the transition dipoles 
of A?2 obtained from HHG data (solid lines) in comparison 
with the theoretical ones (dash lines) in the case of the 800 
nm laser (a) and 1200 nm laser (b). 
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FIG. 6: (Color online) Same as in Fig.[5]but for y-component. 



tirely surprising since one should expect the simple model 
Eq. Q to agree better with Lewenstein model near the 
cut-off, where the contribution from single returns dom- 
inates (Hi. 

In Fig. [7] we show the extracted dipole moments for O2 
molecules using the 1200 nm laser. In this case Xe(5p) 
(with I p = 12.13 eV) is chosen as the reference atom. 
The normalization is done in the similar fashion as for 
N 2 , but here we chose to normalize the y— component of 
the dipole at k 2 = 2.0 for 8 — 45°. We see that the results 
are also comparable with the theoretical data in the range 
of k 2 = [1.0 — 7.5] corresponding to the useful range of 
HHG from H25 to H109. The difference in the cut-off 
positions is due to the different ionization potentials of 
N 2 (15.58 eV) and 2 (12.03 eV). 

It should be noted that for 2 the extracted dipole mo- 
ments for the angles near 0° and 90° do not agree well 
with the theoretical data. It is related to the accuracy 



04 




2 4 6 8 

k 2 (a.u.) 



FIG. 7: (Color online) The transition dipole matrix element of 
O2 extracted from HHG data (solid lines) in comparison with 
the theoretical ones (dash lines) in the case of the 1200 nm 
laser. The upper and lower panels are for y and x components, 
respectively. 



of the small ionization rates and small HHG yields in 
the small angles and near 90°. Note that plane- wave ap- 
proximation is also questionable near these angles, where 
the HOMO of 2 has nodal surfaces. However, it will 
be shown in the next section that these inaccuracies do 
not affect much the wave function obtained by the tomo- 
graphic procedure. 



C. Retrieving wavefunctions by tomographic 
procedures 

We next examine the retrieved wavefunction from the 
extracted "experimental" dipole moment. We emphasize 
that in the following the wave function is understood as 
integrated over z direction (in the molecular frame) as 
given in Eq. (J5J). In Fig. JHJfa.) the contour plot for the 
wavefunction of N 2 obtained from the HHG by the 800 
nm laser is shown. The harmonics used are from H17 to 
H47. In Fig. [5Jb), the wavefunction extracted from the 
HHG by the 1200 nm laser is shown, using harmonics 
orders from H27 to HI 13. The converged wavefunction 
obtained from the tomographic procedure using the theo- 
retical dipole moments in the large range of k 2 = [0—100] 
is shown in Fig. H^c). This is exactly the ground state 
molecular wavefunction used in the calculation of the 
HHG. In the following we will call the HOMOs that are 
obtained from the GAMESS and used to generate HHG 
spectra "exact" wavefunctions. 

By comparing the three figures in Fig. [5J we see that 
the result from the 800nm is not good enough, such that 
even the nodal surfaces are not reproduced. This is due to 
the limited range of momentum space covered by the H17 
to H47 harmonics used in the tomographic imaging. For 
the 1200 nm case, the retrieved wavefunction resembles 
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FIG. 8: (Color online) The contour plot of the wave function 
of N2 obtained from HHG data in the useful range H17 - H47 
of the 800 nm laser (a); from HHG data in the useful range 
H27 - H113 of the 1200 nm laser (b); and from the theoretical 
transition dipole in the range k 2 — [0 — 100] (c). 
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FIG. 9: (Color online) N 2 wave-function "$>(x,y = 0). The 
solid (dashed) curves correspond to the extracted (theo- 
retical) wave-function, respectively. The theoretical wave- 
function are obtained using the tomographic procedure with 
the same range as for the extracted wave-functions. The 
range is k 2 — [1.0 — 7.5] for 1200 nm laser (blue curves) and 
k 2 = [0.9 — 4.2] for 800 nm laser (black curves). For compar- 
ison we also plot the converged wave-function obtained with 
the range k 2 = [0 - 100] (dashed red curve ) . 



the "exact" one quite well, in terms of the orbital shape 
and the nodal surfaces. This clearly illustrates that the 
tomographic procedure will work better if longer wave- 
length lasers are used, such that the plateau region can 
be extended over many harmonics. For more quantita- 
tive comparison, we compare the wavefunctions along the 
internuclear axis derived by using the various approxima- 
tions. In Fig. [9l the "exact" results are compared to the 
wavefunctions obtained from the 800nm and from the 
1200 nm lasers. For each wavelength, we compare the 
wavefunctions derived by using the dipole moments ob- 
tained from the "experimental" HHG data and from the 
theoretical calculations, but within the same range of k 2 . 
Clearly the derived wavefunctions using the same range 
of k 2 are very close to each other. This indicates that 
the dipole moment extracted using the simple three-step 
model, Eq. ([3]), works well at least within the Lewenstein 
model. 

The results in Figs. [5] and [5] indicate that the to- 
mographic imaging method works reasonably well when 
1200 nm laser is used to generate the HHG. We test 
the procedure for two additional internuclear separations, 
R . = 3 a.u. and R — 4 a.u. In Fig. [TU] we compare the 
wavefunctions retrieved from the HHG, using the same 
range of HHG from H27 to H113. Comparing to the 
"exact" wavefunctions, we see that the positions of the 
nodal surfaces, and the symmetry of the wavefunctions, 
are all quite nicely reproduced. Note that in this and the 
following figure for contour plot of the wavefunctions, for 
convenience, we have renormalized the wavefunctions to 
1. It appears that the tomographic procedure fails to cor- 
rectly reproduce the electron density distributions. This 



is especially true for R = 3 and 4 a.u., where it overesti- 
mates the distribution between the two nuclei. We note, 
however, that this is largely due to the limited range of 
harmonics available. The calculations using the theo- 
retical dipoles with the same range of k (or harmonics) 
indeed show the similar trend. 

We further test the tomographic procedure for the O2 
molecule whose HOMO has different symmetry (ir g in- 
stead of a g in N2 case). We used 1200 nm laser to gen- 
erate the HHG spectra. The results for the equilibrium 
distance R = 2.3 and two other internuclear separations 
R = 3 and 4 a.u., are shown in Fig. [TTJ In this case we 
notice that the 7r g orbitals are well reproduced, and the 
internuclear separations, as estimated by the distances 
between the peaks along cc-axis, appeared to be quite ac- 
curate as well. In fact, we have tested the procedures for 
R = 3.5, 4.5 and 5.0 a.u. as well. The deduced internu- 
clear distances from the HHG's are shown in Fig. [121 In 
Table I, the details of the numerical results are shown. 
The input R, the internuclear distance R* read from the 
two peaks of the retrieved wavefunctions are compared. 
When the large range of k 2 (0 to 25 a.u.) in the theo- 
retical dipoles is used in the tomographic procedure, the 
internuclear separations can be extracted to better than 
one percent. Using the smaller range of k 2 of 1.0 to 7.5 
a.u., corresponding to H25 to H109, we obtained R* ac- 
curate to a few percents, except for the R = 2.3 case, 
which is still good to about 10%. 

The relatively good agreement between the retrieved 
and "exact" wave-functions is not entirely surprising. 
This is because Eq. can be considered as an approxi- 
mation to the Lewenstein model used to generate our "ex- 
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FIG. 10: (Color online) Contour plot of the wave-functions 
of N2 with the internuclear distance R = 2.1, 3.0, and 4.0, 
extracted from HHG data (left panels) as compared with the 
exact ones (right panels). 



TABLE I: Internuclear distance of O2 extracted from HHG 
data using the 1200 nm laser. R - internuclear distance as 
input; R* - extracted distance using the tomographic proce- 
dure. The top row corresponds to the case when "experi- 
mental" dipoles are used, the two bottom rows - when the 
theoretical dipoles are used. 
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FIG. 12: Internuclear distance of O2 extracted from HHG 
data using the 1200 nm laser. 
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FIG. 11: (Color online) Same as for Fig. \W\ but for O2 with 
R = 2.3 (top panels), R — 3.0 (middle panels) and R — 4.0 
(bottom panels). 



perimental" data. In particular, this agreement clearly 
cannot be used to justify the plane-wave approximation, 
since this approximation is used in both models. On the 
other hand, one should keep in mind that the Lewenstein 
model has been shown to be able to interpret the main 
features of the measured HHG spectra for atoms and sim- 
ple molecules 0, [2?], HH . We will discuss in details the 
limitations of this approach in the next section. 



IV. INTRINSIC AND PRACTICAL 
LIMITATIONS OF THE TOMOGRAPHIC 
PROCEDURE 

In this section we reexamine the assumptions made in 
the tomographic procedure, following Itatani et al [HI ]. 
The starting point of the method is the Fourier slice 
theorem, Eqs. (6) and (7), which gives a mathematical 
identity between ^(x, y) and the dipole matrix elements 
d(k) in the direction parallel and perpendicular to the 
laser's polarization, with the condition that the transi- 
tion dipole d(k) be given by (k\r\^o(r)), i.e., the final 
continuum states be represented by plane waves. Oper- 
ationally, the dipole matrix elements are to be obtained 
from experiments, for example, using high-order harmon- 
ics generated from aligned molecules, with laser's polar- 
ization direction varying from parallel to perpendicular, 
with respect to the alignment direction of the molecules, 
as carried out by Itatani et al. 

We now examine the assumptions and limitations en- 
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countered in the practical implementation of the tomo- 
graphic imaging method. 

(1) Limitation of the available experimental spectral 
range. 

The Fourier slice theorem requires that the dipole ma- 
trix elements be measured for the whole spectral range, 
while the usable range from the HHG is limited to the 
small plateau region. In Figs. 8-11 we show the effect on 
the retrieved wavefunctions if truncated spectral range is 
used. For the small plateau region available from HHG 
generated by 800 nm lasers, we conclude that accurate 
wavefunctions cannot be obtained, but significant im- 
provement can be made if the HHG generated by 1200 
nm or longer wavelengths are used. We comment that in 
general by introducing a filter function, the accuracy of 
the retrieved HOMO can be improved [l(| [H| . 

(2) Limitation of the approximate three-step model 
for extracting molecular transition dipoles from the mea- 
sured HHG spectra. 

The link between the HHG yields and the "experimen- 
tal" dipole moments of an individual molecule which are 
used in retrieving the molecular wavefunction of the ini- 
tial state is given by Eq. (3) - an equation based on 
the three-step model of the high-order harmonic gener- 
ation. Eq. (3) was not derived from a rigorous HHG 
theory so its validity can only be checked empirically. 
In Eq. (3), propagation effect was included by assuming 
perfect phase matching, an approximation which will be 
addressed later. Given that Eq. (3) is correct, how to 
interpret the derived "experimental" dipole matrix ele- 
ment \d(uj,6)\? In Itatani et al., they interpreted it as 
the dipole matrix element to be used in the Fourier slice 
theorem, Eqs. (6) and (7). This amounts to approximat- 
ing the "exact" continuum states of a molecule by plane 
waves over the whole spectral range, a very poor approx- 
imation well known from photochemistry. 

In this paper, we showed that Eq. (3) is reasonably 
valid for HHG generated from a single molecule if the 
harmonics are calculated using the Lewenstein model, see 
the comparisons in Figs. 5-7. In the Lewenstein model, 
the continuum wavefunction is approximated by plane 
waves, thus the "experimental" dipole matrix element 
from Eq. (3) can be used directly in Eqs. (6) and (7) to 
retrieve the molecular wavefunction. Does Eq. (3) permit 
us to extract accurate dipole matrix elements? We do 
not know yet. Even if it does, it would not help since the 
tomographic procedure employs the condition that the 
continuum states be represented by plane waves. 

A further limit of Eq. (3) is that it gives only the ab- 
solute values of the dipole matrix element. For simple 
linear molecules, the dipole matrix elements are real such 
that one can just change the sign whenever the absolute 
value goes through a near-zero minimum, as adopted by 
Itatani et al. For more complex molecules, the dipole 
moments can become complex and it would be difficult 
to extract the phase information. As for measurement of 
the relative phase of high harmonics, important progress 
has been achieved recently [H, [H, [37j ■ 



(3) Limitations from the macroscopic effects. 

The high-order harmonics measured in the labora- 
tory are from the coherent emission of light from many 
molecules. Eq. (3) was written down under the perfect 
phase matching condition, but how well is it met in actual 
experimental condition, especially from partially aligned 
molecules? We note that in applying Eqs. (6) and (7), 
both parallel and perpendicular components of the HHG 
are needed. This can be done in principle if the molecules 
are fully aligned, i.e., all are pointing in a fixed direction 
in space. Such condition cannot be reached by the pump 
laser pulse where molecules are only partially aligned. 
Thus the angles between the laser's polarization and the 
molecules in the ensemble are not all the same, making 
it very difficult to extract the perpendicular component 
of the dipole moment of the individual molecule from the 
measured HHG signal. Note that for a unaligned ensem- 
ble of molecules, the perpendicular component of HHG 
vanishes. In Itatani et al., the molecules were assumed to 
be fully aligned in the direction of the pump pulse. They 
did not measure the perpendicular component of HHG, 
but the perpendicular component of the dipole moment 
was assumed. It turns out that the form they assumed 
is very close to what we calculated using the Lewenstein 
model, shown in Fig. 2(b). Note that in our present 
model study we did not face the macroscopic issues since 
we derived "experimental" dipole matrix elements from 
HHG due to a single molecule fixed in space only. 



V. CHEMICAL IMAGING WITH HIGH-ORDER 
HARMONIC GENERATION 

A. Alternative methods of chemical imaging 

In Section III.C, with our most optimistic assumptions, 
we have shown that the HOMO wavefunction cannot be 
accurately retrieved based on the tomographic method 
since the spectral range provided by the plateau harmon- 
ics is too limited. On the other hand, the symmetry of 
the HOMO, and the internuclear distances between the 
two atoms, are accurately reproduced from the retrieved 
wavefunction, especially when lasers with longer wave- 
lengths, say, 1200 nm, are used. Such information is com- 
parable to what one can obtain from electron or X-ray 
diffraction experiments. On the other hand, since HHG 
can be generated with femtosecond lasers of duration of 
the order of 10 fs or less, the retrieval of even such limited 
information can be very useful in probing the time depen- 
dence of a chemical transformation at time resolution of 
the order of tens to sub- ten femtoseconds. Thus further 
investigations of possible dynamic chemical imaging with 
HHG by femtosecond lasers are warranted. 

In the future, however, we believe that dynamic chem- 
ical imaging be developed along a different route. As dis- 
cussed in Section IV, the steps involved in obtaining sin- 
gle molecule dipole moments from the macroscopic HHG 
rely on many assumptions which are difficult to remove. 
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At present, we suggest that dynamic chemical imaging 
with HHG be proceeded as follows. We assume that the 
initial configuration of the molecule in the ground state 
is known from the conventional imaging method. Experi- 
mentally the alignment dependent HHG from the ground 
state of such a molecule will first be measured, and the- 
oretical calculations be performed to make sure that the 
experimental data at this initial configuration are well 
described. Under chemical transformation, our goal is to 
locate the positions of the atoms as they change in time, 
by measuring their HHG spectra at different time delays, 
following the initiation of the reaction. Assuming that 
such HHG data are available experimentally, the task of 
dynamic chemical imaging is to find ways to extract the 
intermediate positions of all the atoms during the trans- 
formation, and in particular, to identify the important 
transition states of the reaction. This will be done by 
the iterative method. 

In the iterative method, we will first make a guess on 
the new positions of all the atoms in the molecule. A 
good guess is to follow the reaction coordinates along 
the path where the potential surfaces have the local min- 
imum. Exis ting quantum chemistry codes (see, for ex- 
ample, Ref. [24|, (38|) provide good guidance as a start- 
ing point. For each initial guess the HHG from single 
molecules will be calculated. Then the propagation effect 
and the alignment of molecules are taken into account. 
The resulting macroscopic HHG are then compared to 
the experimental data to find atomic configurations that 
best fit the experimental data. Recall that for each time 
delay, the harmonics can be measured for many different 
angles between the aligning laser and the probe laser. 
Different wavelengths of the lasers can also be used. At 
present, we believe that HHG be measured along the 
probe laser's polarization direction only, since the per- 
pendicular component tends to be severely averaged out 
when molecules are only partially aligned. In principle, 
the phase of each harmonic can also be measured exper- 
imentally [37l ]. For the chemical imaging purpose, how- 
ever, we think it is easier to use different wavelengths 
to generate the larger data set for retrieving the atomic 
configurations instead. We mention that the best fitting 
criterion of the iterative method is still to be developed, 
but at least one can start with the genetic algorithm [40j • 

For the methods suggested above to work, a relatively 
accurate and efficient theoretical model for calculating 
HHG from single molecules is needed. Direct solution of 
the time-dependent Schrodinger equation for molecules 
in a laser field, even within the single active electron 
model is out of questions. On the other hand, a model 
such as the Lewenstein model is very efficient. It is 
straightforward to include the propagation effect and the 
anisotropic distributions of molecules. Thus if the macro- 
scopic HHG's can be calculated efficiently with reason- 
able accuracy, the large set of experimental HHG data 
can then be analyzed to retrieve the atomic configura- 
tions at each time step efficiently. We emphasize that in 
this iterative procedure, the number of parameters that 



are to be determined is small since only the positions of 
the atoms that undergo large changes are to be deter- 
mined. 

We summarize that with this method, there are a num- 
ber of advantages over the tomographic method: (1) No 
need to assume that the continuum states be (incorrectly) 
represented by plane waves; (2) No obvious need to mea- 
sure HHG along both the parallel and the perpendicular 
directions. (3) The method is based on a fitting pro- 
cedure, not a Fourier transformation, such that limited 
spectral range like the plateau harmonics may be ade- 
quate. For the tomographic method the whole spectral 
range is needed. (4) The propagation and the alignment 
effect of the molecules can straightforwardly be incorpo- 
rated in the fitting process. 



B. Outlook of chemical imaging with HHG 

The method we have suggested above for chemical 
imaging formally requires an accurate and efficient the- 
oretical method for calculating the HHG of a complex 
molecule. In a chemical reaction involving "complex" 
molecules, normally only a few atoms make large change 
in positions. Such changes often modify the HOMO's sig- 
nificantly since the symmetry of the molecule is altered, 
resulting in significant change in the alignment depen- 
dence of the HHG. We note that such changes are typ- 
ically accounted for by the Lewenstein model already. 
Thus even before a more accurate theory of HHG be- 
comes available, the procedure outlined in the previous 
subsection can be tested at least qualitatively. 

As an example, consider the isomerization of acetylene 
to vinylidene [4ll, [42. 43]. We have calculated the simpli- 
fied possible reaction path of the isomerization by using 
density functional theory [44j with the hybrid exchange- 
correlation functional of Perdew, Burke and Ernzerhof 
(PBE1PBE) (4f| and Dunning's correlation consistent 
basis set AUG-cc-pVTZ [38( . It is in excellent agreement 
with Ref. [43| . This path involves the well-known transi- 
tion state TS2 [42| . To reach TS2 from acetylene, at least 
about 2 eV is needed. Before reaching the vinylidene 
there is one additional local minimum (LM2) and one 
additional transition state (TS1) in the potential surface. 
The positions of the atoms and the shape of the HOMO 
for each of these special configurations are shown in the 
figure. Note that the binding energy and the HOMO 
for acetylene is very different from the other four. Such 
differences can be easily revealed from the "alignment" 
dependence of HHG. In Fig. 13 we show the HHG cal- 
culated for these different isomers. The alignment angle 
is defined to be between the C-C axis and the laser po- 
larization direction. Except for the linear acetylene, the 
HHG for the other planar molecules should depend on 
the (f) angle which is between the plane of the molecule 
and the plane of the C-C axis and the laser's polariza- 
tion. In Fig. 13, this ^-dependence has been integrated 
over. 
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FIG. 13: (Color online) Predicted angular dependence of the 
21th (red), 25th (green), 29th (blue), 31st (purple) harmonic 
yields in the plateau region for C2H2 isomers: (a) Acetylene; 
(b) TS2; (c) Vinylidene; (d) TS1; and (e) LM2. 
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FIG. 14: (Color online) Predicted angular dependence of the 
21th (red), 25th (green), 29th (blue), 31st (purple) harmonic 
yields in the plateau region for HCN isomers: (a) hydrogen 
cyanide HCN; (b) transition state CHN; and (c) hydrogen 
isocyanide CNH. 



The results in Fig. 13 illustrates the potential and 
the limitations of using the HHG for imaging molecules. 
While it is relatively easy to identify acetylene from 
vinylidene by the alignment dependence of the HHG 
yields, it is difficult to distinguish vinylidene from other 
transient molecules, TS1, TS2 and LM2. Although the 
atomic arrangements for TS2, LM2 and TS1 are quite dif- 
ferent from vinylidene, major features of the HOMO from 
each of these "molecules" are relatively similar. Since the 
HHG, at least within the simplest model, is generated 
from the HOMO only, it is to be understood as a probe 
of the highest occupied molecular orbital, rather than 
a probe of the positions of the atoms in the molecule, 
as from the diffraction method. Thus, in a way, imag- 
ing by HHG is more similar to (e,2e) momentum spec- 
troscopy [4y,|47|. While it may be possible to use HHG to 
distinguish vinylidene from TS2, LM2 and TS1, by care- 
fully analyzing the dependence of HHG on the alignment 
angle 9 and the azimuthal angle (p, or by using lasers with 
different wavelengths, the procedure is likely much more 
demanding. On the other hand, in a typical reaction, a 
vibrational wave packet is created. Since the potential 
surfaces of these four states are quite close to each other, 
unique identification of the atomic configurations of these 
four isomers as stationary states may be an oversimpli- 
fication. Nevertheless, other imaging methods would be 
needed to map out the full evolution from acetylene to 
vinylidene. The transformation away from the acetylene, 
on the other hand, would be very easy to identify using 
the HHG. 

Before closing, we offer one more example - the isomer- 
ization of hydrogen cyanide HCN to hydrogen isocyanide 
HNC [48|. We have calculated the local potential min- 
ima, the shape of the HOMO, the atomic configurations, 
and the contour plot of the potential energy surface cal- 
culated by using DFT [44[ with Becke's three parame- 
ter hybrid functional B3LYP [H and 6-311+G(2df,2pd) 
Gaussian- type basis set [38j |. The latter provides infor- 
mation on the reaction path from the linear HCN to the 
linear HNC, via a transition state of a coplanar molecule 
as sketched. A fair large energy of 2.1 eV is needed to 
reach the transition state. Note that the HOMO's for 
the three "molecules" are rather different, and these dif- 
ferences are easily revealed from the large difference in 
the alignment dependence of the HHG yields shown in 
Fig. 14, calculated again using the Lewenstein model. 
Imaging the transition from HCN to HNC via the tran- 
sition state using the HHG would appear to be much 
simpler in this case. (We comment that in this exam- 
ple, Coulomb explosion technique probably would be the 
easiest since there are only three atoms in the system.) 

In the discussion so far, the molecules have been as- 
sumed to be stationary with its major axis making an 
angle 9 with respect to the laser polarization. Linear 
molecules can be aligned by a short linear-polarized laser, 
or oriented in a combined DC field and a laser [Hcj . For 
more complex molecules, they can be oriented with a 
elliptically polarized laser pulse [39]. Clearly the orien- 
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tation and alignment of the molecules itself is a chal- 
lenging subject on its own. With the rapid progress 
of laser technology, we are optimistic that HHG from 
aligned and/or oriented molecules will become available 
for dynamic imaging in the coming years. 

VI. SUMMARY AND FUTURE DIRECTIONS 

In this paper we examined the theoretical basis of the 
tomographic method used by Itatani et al. [la ]. Due to 
the limited spectral range provided by the high order 
harmonics in the plateau region, we show that intrin- 
sically it is very difficult to reconstruct accurately the 
wavefunction of the highest occupied molecular orbital 
(HOMO) from the measured high order harmonics, but 
the symmetry of the HOMO and the internuclear sep- 
aration can indeed be adequately obtained. The tomo- 
graphic method relies on extracting the dipole moment of 
individual molecules from the HHG measured in a macro- 
scopic medium, and that the continuum wavefunctions in 
the dipole moment of a molecule be represented by plane 
waves. The former posts many practical difficulties owing 
to the macroscopic propagation effect and the imperfect 
alignment of molecules in the laser field. For the latter, 
representing continuum wavefunctions of a molecule by 
plane waves is a drastic simplification. These limitations 
pose severe obstacles in any attempt to improve the to- 
mographic imaging of molecules. 

In spite of these limitations, we showed that useful in- 
formation on the molecule, in particular, the nature of 
the HOMO and the position of the atoms in a transient 
molecule, can be extracted from the measured HHG's. 
In view of the limitations of the tomographic method, we 
have proposed an iterative approach to extract such in- 
formation in the future. The method relics on an efficient 
and relatively accurate theoretical model for calculating 
HHG from molecules. With the availability of such a the- 
ory, the propagation effect and other macroscopic effects 
can be modelled to generate HHG that can be compared 



to the measured data. However, such a simple theory 
for the HHG by molecules is not yet available today. On 
the other hand, the alignment dependence of the HHG 
emission is very sensitive to the nature of the HOMO. 
If the HOMO changes significantly when the molecule 
undergoes transformation, for example, in isomerization, 
then the HHG yield would reveal the changes efficiently. 
Using HHG calculated from the Lewenstein model for 
the isomerization of C2H2 and of HCN, we showed that 
dynamic chemical imaging with HHG is likely possible 
even without detailed quantitative iterative procedures 
outlined in Section V. 

Still many challenges remain for dynamic chemical 
imaging with lasers. Theoretically it is imperative that 
a simple and accurate theory of HHG from molecules 
be developed. Experimentally, techniques for efficiently 
aligning molecules with lasers are needed. Since the 
HHG yield from each molecule is small, the fraction of 
molecules that undergo chemical transformation should 
be large within the target volume. Because HHG can also 
be generated from the unexcited ground-state molecules, 
the intensity of the probe lasers should be selected ap- 
propriately. One advantage of using HHG for chemical 
imaging is that it actually probes the HOMO instead of 
the atomic positions of the whole molecule, thus it is 
much more directly relevant to chemical transformation 
since it is the change of the outermost electron orbital 
that dictates chemical reactions. 
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